Non-equilibrium dynamics of Ising models with decoherence: an exact solution 
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The interplay between interactions and decoherence in many-body systems is of fundamental 
importance in quantum physics: Decoherence can degrade correlations, but can also give rise to a 
variety of rich dynamical and steady- state behaviors. We obtain an exact analytic solution for the 
non-equilibrium dynamics of Ising models with arbitrary interactions and subject to the most general 
form of local Markovian decoherence. Our solution shows that decoherence affects the relaxation 
of observables more than predicted by single-particle considerations. It also reveals a dynamical 
phase transition, specifically a Hopf bifurcation, which is absent at the single-particle level. These 
calculations are applicable to ongoing quantum information and emulation efforts using a variety of 
atomic, molecular, optical, and solid-state systems. 
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Understanding strongly correlated quantum systems in 
the presence of decoherence is a fundamental challenge 
in modern physics. While decoherence generally tends to 
degrade correlations, it is now widely appreciated that it 
can also give rise to many-body physics not possible with 
strictly coherent dynamics p~H6] , and can be used explic- 
itly for the creation of entanglement [TUTU] . Regardless of 
whether one's intention is to minimize or to harness deco- 
herence, determining its effect on interacting many-body 
systems is central to quantum simulation [TT], quantum 
information [12], and quantum metrology [13 j. So far, 
this understanding has been hindered by the computa- 
tional complexity of numerical techniques for open sys- 
tems and the scarcity of exact analytic solutions. Exact 
solutions for dynamics of interacting quantum systems in 
dimensions greater than one are rare even in the absence 
of decoherence, and to our knowledge no such solutions 
have been obtained in the presence of local decoherence. 

The central result of this manuscript is an exact solu- 
tion, Eqs. (9][TT), for the time-dependence of all two-spin 
correlation functions in a system of spins interacting via 
arbitrary Ising couplings, and subject to the most general 
form of local Markovian decoherence allowed by nature 
[T4] . Our solution is applicable to a broad range of impor- 
tant quantum systems, including trapped ions [T5HT7] , 
polar molecules [18] [19] , Rydberg atoms [20] [21] , neutral 
atoms in optical cavities [22] [23], optical lattice clocks 
[24] , superconducting qubits[25 , quantum dots [26 , and 
nitrogen vacancy centers [27]. Here we apply our solution 
to trapped ion experiments because: (1) the relative im- 
portance of decoherence and coherently driven quantum 
correlations is controllable, (2) the tunable long-ranged 
interactions are generically frustrated, making large-scale 
numerical simulations impractical, and (3) these experi- 
ments are the most developed, with non-equilibrium dy- 
namics already being explored [T5] . 



FIG. 1: (Color online) Schematic illustration of the various 
decoherence processes. The figure shows a lattice of spins 
initialized to point in a direction orthogonal to the z-axis. 
Longitudinal (Ti) spin relaxation with rates r u d, Td u (red 
and blue spin, respectively), and dephasing (T2) with rate Y e i 
(green spin) are shown for three different spins. In atomic 
systems, one way this decoherence can arise is due to sponta- 
neous emission from an excited level (top panels) [28] . 



equilibrium dynamics of a long-ranged Ising Hamiltonian 
^ = 77E J ^I- (1) 



i<j 



Here a a (a — x, y, z) are Pauli matrices, M is the to- 
tal number of spins, and Latin letters from the middle 
of the alphabet are site indices. Our results are valid for 
arbitrary , but given the relevance to numerous exper- 
iments we will sometimes consider power law couplings 
Jij = J\ri — where is the position of the i th 

spin in lattice units (( = 3 for polar molecules, ( — 6 
for Rydberg atoms, and < ( < 3 for trapped ions). 
In the presence of local decoherence, the most general 
Markovian dynamics of the system reduced density ma- 
trix obeys a master equation [37] 



Hp=-i (Uesp-pHlg^+Vip). 



(2) 



Outline of the calculation. — We consider far-from- 



The effective Hamiltonian H e ff and dissipator V 
have contributions from all jump operators J G 
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= U-i^J ] J, V(p) = 2^2 JpJ^ (3) |*i f * 2 | *3 f *4 | t 



all J 



all J- 



where = \{cr^ ± ^J)- The jump operators <t~, <r + , 
and a z give rise to spontaneous de-excitation, sponta- 
neous excitation, and elastic dephasing, respectively (see 
Fig. [I]). We refer to the spin-changing processes (a ± ) as 
Raman decoherence, and the spin-preserving processes 
(a z ) as Rayleigh decoherence. In what follows we as- 
sume an initially pure and uncorrelated density operator 

P (0) = |V(0)><v>(0)|, with |v>(0)> = (g^/j^k,). 



Here o\- 



±1 are the eigenvalues of <t|, /j(l) 



cos((9 j /2)e i ^/ 2 , and ^(-1) = sin((9 J /2)e- i ^/ 2 , for ar- 
bitrary Oj and 

Our approach to the solution of Eq. Q for the chosen 
initial conditions and arbitrary {r ud , r du ,r e i} is based 
on the quantum trajectories method [29], in which 1^(0)) 
is time-evolved with the effective Hamiltonian 



^eff — T~L 



rei 
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(4) 



and the dynamics are interspersed with stochastic ap- 
plications of the jump operators. In Eq. Q we have 
defined T r = r ud + r du and 7 = |(r ud -r du ). According 
to the standard prescription [29], a particular trajectory 
consists of a set of jump times {£1,^2, •••}, which are se- 
lected by equating the norm of the wave-function to a 
random number uniformly distributed between and 1. 
Choosing which jump operator to apply at each time re- 
quires calculating all expectation values (J^J). Because 
T-L is Hermitian and commutes with all products 
it has no effect on the selection of the jumps, which can 
therefore be obtained for each spin independently (since 
the anti-Hermitian part of H e ff does not couple differ- 
ent spins). With the jump times and jump operators in 
hand, we define a string of rij (time-labeled) jump oper- 
ators on site j as Qj = Jj{t^) x • • • x J™ 3 (t™ 3 ). The time 
evolution of the wavefunction along a trajectory is then 



\m)=r(e- iH ^i[Q j )\m), 



(5) 



where the time-ordering operator T enforces that the 
jump operators are interspersed in the time evolution ac- 
cording to their time labels. 

The time-ordering of the Rayleigh jumps can be ig- 
nored: Because [<rJ,'Heff] = and [<tJ,<t^] = ±2o^, all 
Rayleigh jumps can be evaluated at t — (their commu- 
tation with Raman jumps affects only the overall sign of 
the wavefunction). To the contrary, the Raman jumps do 
not commute with H e ff, and their time ordering cannot 
be so easily accounted for. However, imagine that spin j 
undergoes a single Raman jump, created by applying &j~ 
at time t. This jump operator not only flips spin j into 



h + (t 3 -t 2 ) + (t-t 4 ) 



FIG. 2: (Color online) A series of Raman flips of the spin 
on site j can be formally accounted for as a magnetic field of 
strength 2J jk /N acting for a period of time rj p - r/ own . In 
the notation defined below, this series of jumps is represented 
by the operator Qj oc (t\)frj (fa)^ (t$)(tj (t±). 



the up position, but also removes all parts of the wave- 
function in which spin j pointed up immediately before 
time t. Hence it is as this spin pointed down before time 
t, and up after time t. Since spin j is always in an eigen- 
state of <jJ, it is a spectator to the Ising dynamics, but it 
does influence the other spins via the Ising coupling J^; 
formally it acts on spin k as an inhomogeneous magnetic 
field of strength 2Jjk/M that pointed down before t and 
up after t. For a spin on site j that undergoes multiple 
Raman processes, the same reasoning allows us to treat 
it as a field of strength 2Jjk/M that acts for a period 
of time Tj = (r/ p - r/ own ), where r ^ down ) i s the total 
amount of time that spin j spends pointing up (down) 
along z (see Fig. [2j. Hence we are free to evaluate all 
of the jump operators at t = to give = ]\- Qjl^), 
thus ignoring the time-ordering in Eq. (J5|, so long as we 
evolve |^) with a modified time-evolution operator 



U = exp 



-it H + J2^j ~ ^Wj 



(6) 



Here rjj = J2 k Jjk^k accounts for the magnetic field of 
all Raman-flipped ions, H is obtained from H by ignor- 
ing the spin-spin couplings to spins that have undergone 
Raman jumps, and 7 accounts for the non- Hermitian part 
of H e ff that is not proportional to the identity operator. 

The expectation value of an arbitrary operator O at 
the end of a particular quantum trajectory is therefore 
given by (6) = $\W OU\$) / $\WU\$) , and formally 
taking the average over all trajectories (denoted by an 

overbar) we have tr[p(9] = (O). Along with Eq. 
these constitute a formal solution for the dynamics of 
any observable. We now proceed to derive closed- form 
expressions for the transverse spin-length and spin-spin 
correlation functions — which have not been derived pre- 
viously even in the absence of decoherence. These are 
the central result of this paper. 

Transverse- spin length and correlation functions. — 
Relaxation of the transverse spin-length in an Ising-type 
spin model is a canonical example of equilibration in a 
closed quantum system [50] F3 2l . In the present model, 
such relaxation occurs due to a combination of the pro- 
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liferation of quantum fluctuations and the equilibration 
with the environment (decoherence). Our theoretical 
treatment allows for both effects to be treated simulta- 
neously, and therefore, in principle, the disentangling of 
these two physically different (but consequentially simi- 
lar) processes. While such relaxation does not directly 
indicate quantum correlations or entanglement, in the 
absence of decoherence it nevertheless is entirely due 
to the buildup of quantum correlations — at the mean- 
field level coherent relaxation is absent. The correla- 
tions that develop during the dynamics can be under- 
stood in more detail by looking at two-spin correlation 
functions — for instance, these furnish a complete descrip- 
tion of spin squeezing [33]. In the absence of deco- 
herence and for all-to-all coupling (£ = 0), it is well 
known [33] that the transverse spin component revives 
at a time r r = Nitfo/{2 J), and that a highly-entangled 
macroscopic-superposition state (MSS) appears at time 
r r /2. This state is characterized by vanishing spin length 
but maximum transverse spin fluctuations, and our solu- 
tion can be used to assess the robustness of such fluctu- 
ations against decoherence. 

To calculate the transverse spin-length along a par- 
ticular trajectory, we assign the discrete- valued vari- 
ables IZj and Fj to each lattice site, which count the 
number of Raman jumps and Rayleigh jumps, respec- 
tively. As we have discussed earlier, all jump opera- 
tors can be applied at t = 0, and therefore specifying 
the {IZj , Tj , Tj} on each site fully determines transverse 
spin-length along that trajectory. The transverse spin 
component in direction ip and on site j is given sim- 
ply in terms of the spin-raising operator on that site, 
(SJ) = cos (p(Sf) + sin<p(S?) = Re[e-^(5+>], and gen- 
eralizing solutions obtained in Ref. [30] [34] we find, 



<*/> 



a j /3 j sin 6^ 
2 5j (2 7 t) 



n 



e 2zJ % Tk gk[2a k t(j ~ iJjk/N)) 



g k (2^ta k ) 

Here g^x) = ^ \fj(a)\ 2 e ax , otj = 5 nj ,o (S being the 
Kronecker-delta symbol), f3j — (—1)^ , and the details of 
the calculation are given in the Supplementary material. 
Defining a function V(7Z, J 7 , r) that determines the prob- 
ability distribution of these variables on a given lattice 
site, we have 



<*7 



E E [■■■ fl[d^nn k ,T k ,T k )(a+). 



(8) 



all 7?. all T 



Equation (8j) constitutes a formal solution for (<5\^), 
and it can always be evaluated efficiently by averaging 
(<7+) over stochastically generated trajectories. However, 
because the noise is uncorrelated from site to site, and 
accordingly the expression inside the product of Eq. ([7]) 
depends only on local stochastic variables, these sums 
and integrals factor into M independent sets (each over 
the three stochastic variables), admitting closed form ex- 



pressions. At this point, to avoid unnecessary complica- 
tions in the ensuing expressions, we will take our initial 
state to point along the x axis (0 = 7r/2, cp = 0), but 
our results easily generalize [35]. In the Supplement we 
explain how to evaluate these sums and integrals, and 
here we simply quote the result. Defining 

$(J,t) = e~ xt cos (t\/s 2 -r) + At sine (t\/ ' s 2 - rj , 

with A = T r /2, s = 2^7 + 2J/JV and r = r ud r du , we find 
1 



-Ft 



Mi 



(9) 



where the total decoherence rate is defined T = ^(T Y + 
r e i) [38 . Similar calculations to those described above 
yield spin-spin correlation functions 



i 



-2Vt 



H QfaJfl + vJuut) (10) 



= \e- vt ^^J jk ,t) n ^J jU tl (11) 



with^(J,t) = e~ xt (is-2j)tsinc(tVs 2 - r) and /i, v = ±. 
These correlation functions, along with similar ones ob- 
tained by interchange of the site indices, completely de- 
termine the spin-spin correlations. Each instance of o x 
or a y in an observable generates an overall multiplica- 
tive factor of e~ Ft ; this is the effect of decoherence at the 
single-particle level. The structure of $(J, t) captures 
the effect of decoherence on the many-body physics, and 
could not have been deduced without our exact treat- 
ment. 

We note that for 7 = the function <!>( J, t) undergoes a 
qualitative transition from oscillatory (s 2 > r) to d amped 



(s 2 < r) behavior when T r = T c r = U/Af [Fig. [3(a) 
Therefore, when there is only one coupling strength J, for 
instance in the case of all-to-all or nearest-neighbor cou- 
plings, dynamics of the transverse spin- length and cor- 
relation functions undergoes the same transition. More 
generally, for nonzero 7 and any correlation functions 
exhibit a series of Hopf bifurcations as T_ L is tuned. We 
note that the factors of e~ rt in Eqs. 



(9- 



11 ) may make 



this transition difficult to observe experimentally, since 
the correlations are rapidly suppressed at the critical T^. 

Application to a trapped ion quantum simulator. — 
Trapped ion systems can simulate the Hamiltonian in Eq. 
([!]), and can accurately measure the decoherence rates 
Teh Tud, and T^ u . We note that both the Born-Markov 
approximation and the assumption of uncorrelated deco- 
herence processes are extremely well justified for trapped 
ion systems [28]. Sample averaged spin-length and spin- 
spin correlation functions are easily measured in these ex- 
periments by looking at the length (and its shot-to-shot 
fluctuations) of various projections of the Bloch vector. 
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FIG. 3: (Color online) (a) Plots of <£( J, t) for 7 = and Y v /Y c v £ {0, 1/4, 1/2, 3/4, 1}, showing a transition from oscillatory to 
damped behavior, (b) Example of how spin squeezing is affected by decoherence (dashed lines are Y = 0.06J, solid lines are 
for r = 0) for long-ranged (red, £ = 0) and short-ranged (blue, C — 3) interactions, (c) Transverse spin relaxation and revivals 
for (" = 0, with parameters corresponding to expected experimental capabilities in [15] (blue solid line). The dotted blue line is 
obtained by treating decoherence at the single-particle level, and underestimates the detrimental effect of Raman decoherence 
by about a factor of 35. Inset: transverse-spin fluctuations peaking at time r r /2. Experimental parameters, exact treatment 
(red solid line); experimental parameters, single-particle treatment (red-dotted line); no decoherence (black dashed line). 



In the trapped ion experiments discussed in Ref. [15], 
when £ = the time scale at which quantum correla- 
tions become important for these observables, r c , scales 
with some power of the ion number: r c ~ A/" 1 / 3 for spin- 
squeezing, ~ A/* 1 / 2 for transverse-spin relaxation, and 
- J\f for the creation of MSS's. Taking J\f = 100 and 
T = 0.06J, as is typical in that experiment, we expect 
the proper incorporation of decoherence to be quantita- 
tively important even for spin-squeezing, despite it being 
a relatively short-time indication of entanglement. 

Equations ( 1Q|1 1 ) allow us to exactly calculate the ef- 
fect of decoherence and the finite range of interactions 
on the maximum spin squeezing achievable in experi- 
ment. Figure |3(b)| shows the expected squeezing and 
anti-squeezing as a function of angle ip in the yz plane for 
C = 0, 3, with T = 0.06 J and T el = 8r ud = 8r du (typical 
experimental numbers in [15 ]). The effects of decoher- 
ence are more pronounced for shorter-range interactions 
due to the longer time scales for maximal squeezing. For 
this calculation the spins are assumed to be initialized 
(prior to the Ising dynamics) in a coherent state point- 
ing along the x-axis, and we define the spin squeezing 
parameter 



(S x ) 



(12) 



with AA = V(^ 2 ) - (A) , S^> = lY,jvf, and a 



sin^)^ + cos( , 0)(j|. In Fig. 3(b) , £ has been calculated 
at the time when maximal squeezing occurs. 

For the current experimental parameters there is es- 
sentially no spin revival, and no indication of a MSS at 
r r /2. However, assuming 97 ions and expected improve- 
ments in the experiment [15] (a roughly 50-fold increase 
in the ratio J/T), in Fig. |3(c)| we show that transverse 
spin revivals begin to appear. In that figure, the dot- 



ted line is obtained by treating the decoherence at the 
single-particle level, which amounts to attaching a de- 
caying exponential e~ vt to the operators a x and & y . We 
note that this treatment would be exact if the decoher- 
ence were only of the Rayleigh type (T r = 0). The solid 
line is the full solution from Eq. ([9|; the large (~ 35- 
fold) discrepancy between the single-particle and exact 
results indicates that our theory will be essential for un- 
derstanding the behavior of this experiment. In the inset 
of Fig. |3(c) we plot the transverse spin fluctuations AS X 
at times near r r /2. The peak in these fluctuations (which 
would achieve unity in the absence of decoherence) is a 
result of strong transverse-spin correlations in an emerg- 
ing MSS, indicating that the expected improvements to 
the experiments will bring within reach the production 
of MSS's of ~ 100 ions in the near future. 

Conclusions. — We have exactly solved the non- 
equilibrium dynamics of arbitrary Ising models in the 
presence of single-particle decoherence. These calcula- 
tions provide a rare glimpse into the exact structure 
of relaxation dynamics in open and strongly interacting 
quantum systems. They will allow for quantitative in- 
sights into experiments studying interacting spin-systems 
that are not perfectly isolated from the environment. In 
particular, they provide a rigorous foundation for bench- 
marking trapped-ion quantum simulators when interac- 
tions are strong and decoherence is non-negligible on ex- 
perimental timescales. 
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I 

Expressions for spin length and correlation functions along a single trajectory 

As in the main text, we write the state of a single spin as ^2 a z fj{&j)\vj) •> where crj is an index that takes on 

the eigenvalues of the operator <rj, fj(l) = cos (Oj/2) e*^'/ 2 , and fj(—l) = sin (6j/2) e~ % Vil 2 . The initial state of the 
entire system is taken to be a direct product of states for each individual spin 

|V(o)> = ®5>te)to> (SI) 

3 a 3 

= AW) x x M^)kfV--^>- 
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As discussed in the main text, evaluating all of the Rayleigh jumps at t = can be accomplished by changing 
ifj —> (fj + ttTj, which rotates the spin on site j by an angle tt if it has undergone an odd number of Rayleigh 
jumps. Raman jumps can be incorporated at t = by setting Oj = 0(tt) if the final Raman jump on site j was 
cr~f(cr~). We therefore define fj{crf) to be the modification of ffcrf) under those transformations. We also note that 
spins having undergone one or more Raman jumps are treated as an effective magnetic field, and not included in the 
spin-spin coupling term of the Hamiltonian. This is accomplished by changing — ^ a^ajJij (with ctj = &fc -,o) m the 
Hamiltonian and including an effective magnetic field. Therefore, we can write the (time-dependent) wavefunction 
\${t)) evolving under U [as defined in Eq. (K3| of the manuscript] as 
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In order to calculate the transverse spin-length we will first calculate = (^{t)\a^\^{t)) / {^(t)\^{t)) , and at the 

end obtain (S x ) = Re (£,-(5+)). 

Let's imagine, in particular, calculating |<jf \ip(t)} (there is nothing special about the first spin, this just makes 
the notation in what follows less confusing). Because the wavefunction enters twice, this would involve two sums like 
the one in Eq. (S2), over <rj and erf, but very few terms survive: We need er\ = —1, erf = 1, and for all j / 1 we 

must have cr| = erf , so the matrix element is given by 
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If aj =0, the j th spin always has a well defined value of cr| and the choice to include the term 7<j| in the exponentiated 
sum (or not) affects only the overall normalization of the wavefunction. By multiplying 7 by ctj (in the exponentiated 
sum), we have chosen to not include the term 70-J, and this is properly accounted for when normalizing this expectation 
value below. In order to obtain we must divide by the (non-conserved) normalization of the wavefunction 

(tp(t)\tp(t)). Defining gj(x) = £ a \fj( CT )\ 2e ~ ax ( as m the manuscript), we obtain 
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This expression can be simplified by making the following set of observations: (1) /* (l)/i(— 1) = ai/^ie^ 1 sin(#i)/2 
[as in the text, /3j = (—1)^], (2) the ol\ in the exponent is irrelevant, because if it takes the value the entire 
expression for(^(t)|<7^| ?/>(£)) vanishes, and (3) the summand factorizes into a product where each term contains only 
local (i.e. on a single site) variables, and hence the sum of products can be exchanged for a product of sums. Taking 
(1-3) into account we obtain 
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Therefore, we can write 
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The calculation of correlation functions follows from extremely similar considerations. For instance, let's consider 
Qff = (crjcr~^). In this case, the operators in the expectation value only connect two states if —erf = cr| = 1, 

cr^f = — erf = 1, and erf = erf whenever I £ {j,k}. Therefore, much as before we have 
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Computing correlation functions involving a single a z , such as C/j^ = {^cr k ) can be achieved by inserting <rj into 
the sum in Eq. (S4), yielding 
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Assuming one or more Raman flip occurred, the variable ftj takes on the values ±1 if the final Raman jump is a - 



Analytic evaluation of stochastic averaging of trajectories 



At this point in the calculation, for clarity of presentation, we set cpj = and 3 = tt/2 (for all j), so all spins point 
along the x axis at t = 0. Defining V(H, r) to be the probability distribution of the variables 7£, J 7 , and r on a 
single site (it is the same on every site), the trajectory averaged expectation value is given by 
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To begin, we note that the probability distribution can be decomposed as V(TZ, J 7 , r) = Ve\(F)V v (lZ,, t), which is valid 
because the probability of Rayleigh jump is independent of whether a Raman jump has occurred (and vice versa). 
The occurrence of random processes follows a Poissonian distribution, so T^F) = e~ Telt / A (Y e \t/ 4)^ / and we have 
also calculated V r . The result depends on whether 1Z is even or odd (the proof is simple but requires some careful 
reasoning, and we do not give it here). We parameterize the 7^-odd solution by fi = (7Z — l)/2 (which will run over 
all non- negative integers), and we parameterize the 7^-even solutions by /i = (JZ — 2)/2 (once again running \i over 
all nonnegative integers), and obtain 
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By evaluating the sum 
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we obtain 
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Because a 3 gives 1 if there have not been any Raman flips at site j and otherwise, the term in the first square bracket 
is just the probability that there has been no Raman flip on site j, which is given by e _Fr ^ 2 cosh(27^) [this comes 
from evolving the wavefunction of a single spin pointing along x with the effective Hamiltonian Eq. Q]. Therefore 
we have 
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Defining s = 2ij + 2 J/A/", the trick is now to evaluate the quantity 
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The second equality follows from pulling off the 1Z = term in the sum, which represents the probability of having 
no Raman flip (and so we can set r = and a = 1 in this term). The second term represents the probability for any 
finite number of Raman flips, and hence we must keep r arbitrary but can set a = 0. The integral over r can be 
evaluated by using the identity 
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where j is a spherical Bessel function. The remaining sum over 1Z can be recognized as a generating function for the 
spherical Bessel functions (or a derivative thereof). Defining parameters A = T r /2 and r = T^T^u, and functions 
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we obtain 



J, t) = e~ M cos(st) + \te~ xt F(st, rt) + ste~ xt G(st, rt) 
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We can now write out the exact solution 
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Because Eqs. ( S7|S8 ) have such a similar structure, the stochastic averaging of correlation functions is almost identical, 
and leads to the similar expressions given in the main text [ Eqs. ([TQ|11 )]. 



